library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(base)
library(xts)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

Data Preprocessing

dataPath <- "/Users/bianca/Documents/U of Chicago/Classes/Time Series"
contracts_volume <- read.csv(paste(dataPath,'Contracts_Volume.csv',sep = '/'), header=TRUE)
contracts_class <- read.csv(paste(dataPath,'Contracts_Classification.csv',sep = '/'), header=TRUE)

Computing Floor Volume

contracts_volume$Total.Volume <- as.numeric(contracts_volume$Total.Volume)
contracts_volume$Electronic.Volume <- as.numeric(gsub(",", "", contracts_volume$Electronic.Volume))

contracts_volume$Floor.Volume <- contracts_volume$Total.Volume - contracts_volume$Electronic.Volume

Join Contracts Class and Volume

joined_data <- contracts_volume %>%
  left_join(contracts_class, by = c("Commodity.Indicator" = "Commodity.Code"))

Sort Volume Data Based on Division

cme_data <- joined_data
cme_data$Division <- "CME"
cme_data$Date <- as.Date(cme_data$Date, format= "%m/%d/%Y")

imm_data <- filter(joined_data, Division == "IMM")
imm_data$Date <- as.Date(imm_data$Date, format= "%m/%d/%Y")

iom_data <- filter(joined_data, Division == "IOM")
iom_data$Date <- as.Date(iom_data$Date, format= "%m/%d/%Y")

Aggregate Volume by Commodity Code and Month

cme_agg <- cme_data %>%
  group_by(Commodity.Indicator, Date) %>%
  summarise(Total_Electronic.Volume = sum(Electronic.Volume),
            Total_Total.Volume = sum(Total.Volume),
            Total_Floor.Volume = sum(Floor.Volume))

imm_agg <- imm_data %>%
  group_by(Commodity.Indicator, Date) %>%
  summarise(Total_Electronic.Volume = sum(Electronic.Volume),
            Total_Total.Volume = sum(Total.Volume),
            Total_Floor.Volume = sum(Floor.Volume))

iom_agg <- iom_data %>%
  group_by(Commodity.Indicator, Date) %>%
  summarise(Total_Electronic.Volume = sum(Electronic.Volume),
            Total_Total.Volume = sum(Total.Volume),
            Total_Floor.Volume = sum(Floor.Volume))

print(cme_agg)
## # A tibble: 21,048 × 5
## # Groups:   Commodity.Indicator [559]
##    Commodity.Indicator Date       Total_Electronic.Volume Total_Total.Volume
##    <chr>               <date>                       <dbl>              <dbl>
##  1 16E                 2011-02-01                      11                 11
##  2 16E                 2011-03-01                      24                 24
##  3 16E                 2011-04-01                      85                 85
##  4 16E                 2011-05-01                      63                 63
##  5 16E                 2011-06-01                      70                 70
##  6 16E                 2011-07-01                      44                 44
##  7 16E                 2011-08-01                       8                  8
##  8 16E                 2011-09-01                      12                 12
##  9 16E                 2011-10-01                       5                  5
## 10 1A                  2007-05-01                       0                116
## # ℹ 21,038 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>
print(imm_agg)
## # A tibble: 5,723 × 5
## # Groups:   Commodity.Indicator [78]
##    Commodity.Indicator Date       Total_Electronic.Volume Total_Total.Volume
##    <chr>               <date>                       <dbl>              <dbl>
##  1 16E                 2011-02-01                      11                 11
##  2 16E                 2011-03-01                      24                 24
##  3 16E                 2011-04-01                      85                 85
##  4 16E                 2011-05-01                      63                 63
##  5 16E                 2011-06-01                      70                 70
##  6 16E                 2011-07-01                      44                 44
##  7 16E                 2011-08-01                       8                  8
##  8 16E                 2011-09-01                      12                 12
##  9 16E                 2011-10-01                       5                  5
## 10 36E                 2011-03-01                      27                 27
## # ℹ 5,713 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>
print(iom_agg)
## # A tibble: 16,399 × 5
## # Groups:   Commodity.Indicator [475]
##    Commodity.Indicator Date       Total_Electronic.Volume Total_Total.Volume
##    <chr>               <date>                       <dbl>              <dbl>
##  1 1A                  2007-05-01                       0                116
##  2 1A                  2008-01-01                       0                 25
##  3 1A                  2009-04-01                       0                 14
##  4 1A                  2011-03-01                      76                426
##  5 1A                  2011-06-01                      26                 26
##  6 1A                  2011-07-01                       0                220
##  7 1A                  2011-08-01                      11                 11
##  8 1A                  2011-09-01                       4                  4
##  9 1A                  2011-11-01                     146                146
## 10 1A                  2011-12-01                      25                 28
## # ℹ 16,389 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>

Aggregate by Month only

cme_agg2 <- cme_data %>%
  group_by(Date) %>%
  summarise(Total_Electronic.Volume = sum(Electronic.Volume),
            Total_Total.Volume = sum(Total.Volume),
            Total_Floor.Volume = sum(Floor.Volume))

imm_agg2 <- imm_data %>%
  group_by(Date) %>%
  summarise(Total_Electronic.Volume = sum(Electronic.Volume),
            Total_Total.Volume = sum(Total.Volume),
            Total_Floor.Volume = sum(Floor.Volume))

iom_agg2 <- iom_data %>%
  group_by(Date) %>%
  summarise(Total_Electronic.Volume = sum(Electronic.Volume),
            Total_Total.Volume = sum(Total.Volume),
            Total_Floor.Volume = sum(Floor.Volume))

Loading Price Data

imputed_cme <- read.csv(paste(dataPath,'imputed_cme.csv',sep = '/'), header=TRUE)
names(imputed_cme) <- c("Date", "Price")
imputed_cme$Date <- as.Date(imputed_cme$Date)
  
imputed_imm <- read.csv(paste(dataPath,'imputed_imm.csv',sep = '/'), header=TRUE)
names(imputed_imm) <- c("Date", "Price")
imputed_imm$Date <- as.Date(imputed_imm$Date)

imputed_iom <- read.csv(paste(dataPath,'imputed_iom.csv',sep = '/'), header=TRUE)
names(imputed_iom) <- c("Date", "Price")
imputed_iom$Date <- as.Date(imputed_iom$Date)

EDA

cme_summary <- left_join(imputed_cme, cme_agg2, by = "Date")
imm_summary <- left_join(imputed_imm, imm_agg2, by = "Date")
iom_summary <- left_join(imputed_iom, iom_agg2, by = "Date")

cme_numeric_data <- cme_summary %>%
  select_if(is.numeric)
cme_cor_matrix <- cor(cme_numeric_data, use = "complete.obs", method = "pearson")
cme_cor <- as.data.frame(cme_cor_matrix)


imm_numeric_data <- imm_summary %>%
  select_if(is.numeric)
imm_cor_matrix <- cor(imm_numeric_data, use = "complete.obs", method = "pearson")
imm_cor <- as.data.frame(imm_cor_matrix)


iom_numeric_data <- iom_summary %>%
  select_if(is.numeric)
iom_cor_matrix <- cor(iom_numeric_data, use = "complete.obs", method = "pearson")
iom_cor <- as.data.frame(iom_cor_matrix)

print(cme_cor)
##                               Price Total_Electronic.Volume Total_Total.Volume
## Price                    1.00000000               0.6530748          0.6802748
## Total_Electronic.Volume  0.65307484               1.0000000          0.9736634
## Total_Total.Volume       0.68027480               0.9736634          1.0000000
## Total_Floor.Volume      -0.06265014              -0.3726483         -0.1512653
##                         Total_Floor.Volume
## Price                          -0.06265014
## Total_Electronic.Volume        -0.37264828
## Total_Total.Volume             -0.15126533
## Total_Floor.Volume              1.00000000
print(imm_cor)
##                             Price Total_Electronic.Volume Total_Total.Volume
## Price                   1.0000000               0.3287382          0.4764459
## Total_Electronic.Volume 0.3287382               1.0000000          0.9465020
## Total_Total.Volume      0.4764459               0.9465020          1.0000000
## Total_Floor.Volume      0.2972396              -0.4706299         -0.1607259
##                         Total_Floor.Volume
## Price                            0.2972396
## Total_Electronic.Volume         -0.4706299
## Total_Total.Volume              -0.1607259
## Total_Floor.Volume               1.0000000
print(iom_cor)
##                               Price Total_Electronic.Volume Total_Total.Volume
## Price                    1.00000000             -0.02515719         0.05596529
## Total_Electronic.Volume -0.02515719              1.00000000         0.98393410
## Total_Total.Volume       0.05596529              0.98393410         1.00000000
## Total_Floor.Volume       0.43943993             -0.29485138        -0.11951925
##                         Total_Floor.Volume
## Price                            0.4394399
## Total_Electronic.Volume         -0.2948514
## Total_Total.Volume              -0.1195193
## Total_Floor.Volume               1.0000000

It looks like Floor Volume is not a good indicator for CME class seats prices. Both Electronic Volume and Total Volume are not good indicators for IOM class seats prices. We can also see that Electronic Volume and Total Volume are highly correlated in all seat classes, so only one of those will be used as predictor.

Task A: Algorithms

Split train and test

train_data_cme <- subset(cme_summary, Date < "2013-01-01")
test_data_cme <- subset(cme_summary, Date >= "2013-01-01")

train_data_imm <- subset(imm_summary, Date < "2013-01-01")
test_data_imm <- subset(imm_summary, Date >= "2013-01-01")

train_data_iom <- subset(iom_summary, Date < "2013-01-01")
test_data_iom <- subset(iom_summary, Date >= "2013-01-01")

price_cme <- ts(train_data_cme$Price, frequency = 12)
price_imm <- ts(train_data_imm$Price, frequency = 12)
price_iom <- ts(train_data_iom$Price, frequency = 12)

Linear Regression

lm_cme <- lm(Price ~ Total_Total.Volume, data = train_data_cme)
lm_forecast_cme <- predict(lm_cme, newdata = test_data_cme, n.ahead = 12)
summary(lm_cme)
## 
## Call:
## lm(formula = Price ~ Total_Total.Volume, data = train_data_cme)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -419623 -112256   -7097   82349  791934 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.166e+05  4.127e+04   2.826   0.0054 ** 
## Total_Total.Volume 2.945e-03  2.371e-04  12.420   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 197600 on 142 degrees of freedom
## Multiple R-squared:  0.5207, Adjusted R-squared:  0.5173 
## F-statistic: 154.3 on 1 and 142 DF,  p-value: < 2.2e-16
lm_imm <- lm(Price ~ Total_Total.Volume, data = train_data_imm)
lm_forecast_imm <- predict(lm_imm, newdata = test_data_imm, n.ahead = 12)

lm_iom <- lm(Price ~ Total_Floor.Volume, data = train_data_iom)
lm_forecast_iom <- predict(lm_iom, newdata = test_data_iom, n.ahead = 12)

For CME linear regression, Electronic Volume and Total Volume are correlated so only Total Volume is used. Floor Volume doesn’t have high correltion with price so it’s not used as predictor. For IMM linear regression, Electronic Volume and Floor Volume is picked instead so there’s no overlap between Total Volume and Floor Volume as predictors. For IOM, only Floor Volume is highly correlated with price so it’s the only predictor.

Linear regression with ARMA errors (use arima with xreg)

arma_model_cme <- auto.arima(price_cme, xreg = train_data_cme$Total_Total.Volume)
arma_forecast_cme <- forecast(arma_model_cme, xreg = test_data_cme$Total_Total.Volume, h=12)

arma_model_imm <- auto.arima(price_imm, xreg = train_data_imm$Total_Total.Volume)
arma_forecast_imm <- forecast(arma_model_imm, xreg = test_data_imm$Total_Total.Volume, h=12)

arma_model_iom <- auto.arima(price_iom, xreg = train_data_iom$Total_Total.Volume)
arma_forecast_iom <- forecast(arma_model_iom, xreg = test_data_iom$Total_Total.Volume, h=12)

Holt-Winters

library(forecast)

hw_model_cme <- hw(price_cme, xreg = train_data_cme$Total_Total.Volume)
hw_forecast_cme <- forecast(hw_model_cme, xreg = test_data_cme$Total_Total.Volume, h = 12)

hw_model_imm <- hw(price_imm,xreg = train_data_imm$Total_Total.Volume)
hw_forecast_imm <- forecast(hw_model_imm, xreg = test_data_imm$Total_Total.Volume, h = 12)

hw_model_iom <- hw(price_iom, xreg = train_data_iom$Total_Total.Volume)
hw_forecast_iom <- forecast(hw_model_iom, xreg = test_data_iom$Total_Total.Volume, h = 12)

ARIMA

arima_cme <- auto.arima(price_cme, seasonal = FALSE)
arima_forecast_cme <- forecast(arima_cme, h = 12)

arima_imm <- auto.arima(price_imm, seasonal = FALSE)
arima_forecast_imm <- forecast(arima_imm, h = 12)

arima_iom <- auto.arima(price_iom, seasonal = FALSE)
arima_forecast_iom <- forecast(arima_iom, h = 12)

Seasonal ARIMA

# price_cme, price_imm, price_iom already has frequency of 12
sarima_cme <- auto.arima(price_cme, seasonal = TRUE)
sarima_forecast_cme <- forecast(sarima_cme, h = 12)

sarima_imm <- auto.arima(price_imm, seasonal = TRUE)
sarima_forecast_imm <- forecast(sarima_imm, h = 12)

sarima_iom <- auto.arima(price_iom, seasonal = TRUE)
sarima_forecast_iom <- forecast(sarima_iom, h = 12)

Fractional ARIMA

# Looking at ACF
acf(price_cme)

acf(price_imm)

acf(price_iom)

The autocorrelations show trend, so arfima is applicable.

library(fracdiff)

d_cme <- fracdiff(price_cme)
st_cme <- diffseries(price_cme, d_cme$d)
arfima_cme <- auto.arima(st_cme, xreg = train_data_cme$Total_Total.Volume)
arfima_forecast_cme <- forecast(arfima_cme, xreg = test_data_cme$Total_Total.Volume, h = 12)

d_imm <- fracdiff(price_imm)
st_imm <- diffseries(price_imm, d_imm$d)
arfima_imm <- auto.arima(st_imm, xreg = train_data_imm$Total_Total.Volume)
arfima_forecast_imm <- forecast(arfima_imm, xreg = test_data_imm$Total_Total.Volume, h = 12)

d_iom <- fracdiff(price_iom)
st_iom <- diffseries(price_iom, d_iom$d)
arfima_iom <- auto.arima(st_iom, xreg = train_data_iom$Total_Total.Volume)
arfima_forecast_iom <- forecast(arfima_iom, xreg = test_data_iom$Total_Total.Volume, h = 12)

ARMA and GARCH

#install.packages("fGarch")
library(fGarch)
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")
num_train_data_cme <- scale(select_if(train_data_cme, is.numeric))
num_train_data_imm <- scale(select_if(train_data_imm, is.numeric))
num_train_data_iom <- scale(select_if(train_data_iom, is.numeric))

num_test_data_cme <- scale(select_if(test_data_cme, is.numeric))
num_test_data_imm <- scale(select_if(test_data_imm, is.numeric))
num_test_data_iom <- scale(select_if(test_data_iom, is.numeric))

garch_model_cme <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_cme, trace = FALSE)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
garch_forecast_cme <- predict(garch_model_cme, data = num_test_data_cme, n.ahead = 12)

garch_model_imm <- garchFit(Price  ~ arma(1,1) + garch(1,1), data = num_train_data_imm, trace = FALSE)
garch_forecast_imm <- predict(garch_model_imm, data = num_test_data_imm, n.ahead = 12)

garch_model_iom <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_iom, trace = FALSE)
garch_forecast_iom <- predict(garch_model_iom, data = num_test_data_iom, n.ahead = 12)

Compare results sMAPE across models

cme_rmse <- data.frame(Commodity = "CME",
                  LM = accuracy(lm_forecast_cme, test_data_cme$Price)[[2]],
                  ARMA = accuracy(arma_forecast_cme, test_data_cme$Price)[[2]],
                  HW = accuracy(hw_forecast_cme, test_data_cme$Price)[[2]],
                  ARIMA = accuracy(arima_forecast_cme, test_data_cme$Price)[[2]],
                  SARIMA = accuracy(sarima_forecast_cme, test_data_cme$Price)[[2]],
                  ARFIMA = accuracy(arfima_forecast_cme, test_data_cme$Price)[[2]],
                  GARCH = accuracy(garch_forecast_cme$meanForecast, num_test_data_cme[,1])[[2]])

imm_rmse <- data.frame(Commodity = "IMM",
                  LM = accuracy(lm_forecast_imm, test_data_imm$Price)[[2]],
                  ARMA =  accuracy(arma_forecast_imm, test_data_imm$Price)[[2]],
                  HW = accuracy(hw_forecast_imm, test_data_imm$Price)[[2]],
                  ARIMA = accuracy(arima_forecast_imm, test_data_imm$Price)[[2]],
                  SARIMA = accuracy(sarima_forecast_imm, test_data_imm$Price)[[2]],
                  ARFIMA = accuracy(arfima_forecast_imm, test_data_imm$Price)[[2]],
                  GARCH = accuracy(garch_forecast_imm$meanForecast, num_test_data_imm[,1])[[2]])

iom_rmse <- data.frame(Commodity = "IOM",
                  LM = accuracy(lm_forecast_iom, test_data_iom$Price)[[2]],
                  ARMA = accuracy(arma_forecast_iom, test_data_iom$Price)[[2]],
                  HW = accuracy(hw_forecast_iom, test_data_iom$Price)[[2]],
                  ARIMA = accuracy(arima_forecast_iom, test_data_iom$Price)[[2]],
                  SARIMA = accuracy(sarima_forecast_iom, test_data_iom$Price)[[2]],
                  ARFIMA = accuracy(arfima_forecast_iom, test_data_iom$Price)[[2]],
                  GARCH = accuracy(garch_forecast_iom$meanForecast, num_test_data_iom[,1])[[2]])

combined_rmse <- rbind(cme_rmse, imm_rmse, iom_rmse)

combined_rmse$Best_Model <- apply(combined_rmse[, -1], 1, function(x) {
  colnames(combined_rmse[, -1])[which.min(abs(x))]
})


print(combined_rmse)
##   Commodity       LM         ARMA         HW     ARIMA    SARIMA    ARFIMA
## 1       CME 271928.7 -241806.2206  -192.5793 31291.667 31291.667 512777.13
## 2       IMM 242626.0   25469.1639 33720.7608 34281.914 34281.914 268904.23
## 3       IOM 116316.6    -900.4912 14852.5595 -1071.768 -4836.629  89559.59
##      GARCH Best_Model
## 1 1.064806      GARCH
## 2 1.242325      GARCH
## 3 1.606774      GARCH

Best Model: - GARCH for all

Task B: BSTS

Load data

sts <- read.csv(paste(dataPath,'sts_data_2023.csv',sep = '/'), header=TRUE)
sts <- na.omit(sts)

sts_1990 <- sts %>%
  filter(as.numeric(format(as.Date(Date, "%m/%d/%y"), "%y")) > 23) %>%
  mutate(Date = format(as.Date(Date, "%m/%d/%y"), "19%y-%m-%d"))

sts_2000 <- sts %>%
  filter(as.numeric(format(as.Date(Date, "%m/%d/%y"), "%y")) <= 23) %>%
  mutate(Date = format(as.Date(Date, "%m/%d/%y"), "20%y-%m-%d"))

sts <- rbind(sts_1990, sts_2000)

sts <- xts(sts[,-1], order.by = as.Date(sts$Date))

#install.packages("bsts")
library(bsts)
## Loading required package: BoomSpikeSlab
## Loading required package: Boom
## 
## Attaching package: 'Boom'
## The following object is masked from 'package:stats':
## 
##     rWishart
## 
## Attaching package: 'BoomSpikeSlab'
## The following object is masked from 'package:stats':
## 
##     knots
## 
## Attaching package: 'bsts'
## The following object is masked from 'package:BoomSpikeSlab':
## 
##     SuggestBurn
data_fit1 <- window(sts, end = "1979-06-30")
data_fit2 <- window(sts,end = "2004-06-30")
data_fit3 <- window(sts, end = "2022-12-31")

new_data1 <- head(window(sts, start = "1979-09-30"), 5)
new_data1_CL <- new_data1[, -3] 
new_data1_GC <- new_data1[, -1] 

new_data2 <- head(window(sts, start = "2004-09-30"), 5)
new_data2_CL <- new_data1[, -3] 
new_data2_GC <- new_data1[, -1] 

new_data3 <- head(window(sts, start = "2023-03-31"), 1)
new_data_3CL <- new_data1[, -3] 
new_data3_GC <- new_data1[, -1] 

Fit Models for CL

Dynamic regression of other 3 variables (lags)

# FIT 1
ss <- AddLocalLevel(list(), data_fit1$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit1$Mean_CL ~ data_fit1$Mean_GC + data_fit1$Mean_CPI + data_fit1$Mean_SPX)

model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit1, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:43 2024 =-=-=-=-=
plot(model_dynamic)

pred_CL_dynamic_1 <- predict.bsts(model_dynamic, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_dynamic_1)

# FIT 2
ss <- AddLocalLevel(list(), data_fit2$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit2$Mean_CL ~ data_fit2$Mean_GC + data_fit2$Mean_CPI + data_fit2$Mean_SPX)

model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit2, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:44 2024 =-=-=-=-=
plot(model_dynamic)

pred_CL_dynamic_2 <- predict.bsts(model_dynamic, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_dynamic_2)

# FIT 3
ss <- AddLocalLevel(list(), data_fit3$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit3$Mean_CL ~ data_fit3$Mean_GC + data_fit3$Mean_CPI + data_fit3$Mean_SPX)

model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit3, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_dynamic)

pred_CL_dynamic_3 <- predict.bsts(model_dynamic, newdata = new_data3, horizon = 1, quantiles = c(.025, .975))
plot(pred_CL_dynamic_3)

Different trend component for each model

# FIT 1
ss <- AddLocalLinearTrend(list(), data_fit1$Mean_CL)

model_trend <- bsts(data_fit1$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_trend)

pred_CL_trend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_1)

# FIT 2
ss <- AddLocalLinearTrend(list(), data_fit2$Mean_CL)

model_trend <- bsts(data_fit2$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_trend)

pred_CL_trend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_2)

# FIT 3
ss <- AddLocalLinearTrend(list(), data_fit3$Mean_CL)

model_trend <- bsts(data_fit3$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_trend)

pred_CL_trend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_3)

State-specification with semi-local linear trend

# FIT 1
ss <- AddSemilocalLinearTrend(list(), data_fit1$Mean_CL)

model_trend <- bsts(data_fit1$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=
plot(model_trend)

pred_CL_semitrend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_1)

# FIT 2
ss <- AddSemilocalLinearTrend(list(), data_fit2$Mean_CL)

model_trend <- bsts(data_fit2$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=
plot(model_trend)

pred_CL_semitrend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_2)

# FIT 3
ss <- AddSemilocalLinearTrend(list(), data_fit3$Mean_CL)

model_trend <- bsts(data_fit3$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=
plot(model_trend)

pred_CL_semitrend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_3)

Fit Models for GC

Dynamic regression of other 3 variables (lags)

# FIT 1
ss <- AddLocalLevel(list(), data_fit1$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit1$Mean_GC ~ data_fit1$Mean_CL + data_fit1$Mean_CPI + data_fit1$Mean_SPX)

model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit1, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:47 2024 =-=-=-=-=
plot(model_dynamic)

pred_GC_dynamic_1 <- predict.bsts(model_dynamic, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_dynamic_1)

# FIT 2
ss <- AddLocalLevel(list(), data_fit2$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit2$Mean_GC ~ data_fit2$Mean_CL + data_fit2$Mean_CPI + data_fit2$Mean_SPX)

model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit2, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:47 2024 =-=-=-=-=
plot(model_dynamic)

pred_GC_dynamic_2 <- predict.bsts(model_dynamic, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_dynamic_2)

# FIT 3
ss <- AddLocalLevel(list(), data_fit3$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit3$Mean_GC ~ data_fit3$Mean_CL + data_fit3$Mean_CPI + data_fit3$Mean_SPX)

model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit3, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=
plot(model_dynamic)

pred_GC_dynamic_3 <- predict.bsts(model_dynamic, newdata = new_data3, horizon = 1, quantiles = c(.025, .975))
plot(pred_GC_dynamic_3)

Different trend component for each model

# FIT 1
ss <- AddLocalLinearTrend(list(), data_fit1$Mean_GC)

model_trend <- bsts(data_fit1$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=
plot(model_trend)

pred_GC_trend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_1)

# FIT 2
ss <- AddLocalLinearTrend(list(), data_fit2$Mean_GC)

model_trend <- bsts(data_fit2$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=
plot(model_trend)

pred_GC_trend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_2)

# FIT 3
ss <- AddLocalLinearTrend(list(), data_fit3$Mean_GC)

model_trend <- bsts(data_fit3$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)

pred_GC_trend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_3)

State-specification with semi-local linear trend

# FIT 1
ss <- AddSemilocalLinearTrend(list(), data_fit1$Mean_GC)

model_trend <- bsts(data_fit1$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)

pred_GC_semitrend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_1)

# FIT 2
ss <- AddSemilocalLinearTrend(list(), data_fit2$Mean_GC)

model_trend <- bsts(data_fit2$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)

pred_GC_semitrend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_2)

# FIT 3
ss <- AddSemilocalLinearTrend(list(), data_fit3$Mean_GC)

model_trend <- bsts(data_fit3$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)

pred_GC_semitrend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_3)

Testing RMSE

CL_fit1_rmse <- data.frame(Commodity = "CL Fit 1",
                  DynamicRegression = accuracy(pred_CL_dynamic_1$mean, new_data1$Mean_CL)[[2]],
                  LinearTrend = accuracy(pred_CL_trend_1$mean, new_data1$Mean_CL)[[2]],
                  SemiLocalTrend = accuracy(pred_CL_semitrend_1$mean, new_data1$Mean_CL)[[2]])
                  
CL_fit2_rmse <- data.frame(Commodity = "CL Fit 2",
                  DynamicRegression = accuracy(pred_CL_dynamic_2$mean, new_data2$Mean_CL)[[2]],
                  LinearTrend = accuracy(pred_CL_trend_2$mean, new_data2$Mean_CL)[[2]],
                  SemiLocalTrend = accuracy(pred_CL_semitrend_1$mean, new_data2$Mean_CL)[[2]])

CL_fit3_rmse <- data.frame(Commodity = "CL Fit 3",
                  DynamicRegression = accuracy(pred_CL_dynamic_3$mean, new_data3$Mean_CL)[[2]],
                  LinearTrend = accuracy(pred_CL_trend_3$mean, new_data3$Mean_CL)[[2]],
                  SemiLocalTrend = accuracy(pred_CL_semitrend_3$mean, new_data3$Mean_CL)[[2]])

GC_fit1_rmse <- data.frame(Commodity = "GC Fit 1",
                  DynamicRegression = accuracy(pred_GC_dynamic_1$mean, new_data1$Mean_GC)[[2]],
                  LinearTrend = accuracy(pred_GC_trend_1$mean, new_data1$Mean_GC)[[2]],
                  SemiLocalTrend = accuracy(pred_GC_semitrend_1$mean, new_data1$Mean_GC)[[2]])
                  
GC_fit2_rmse <- data.frame(Commodity = "GC Fit 2",
                  DynamicRegression = accuracy(pred_GC_dynamic_2$mean, new_data2$Mean_GC)[[2]],
                  LinearTrend = accuracy(pred_GC_trend_2$mean, new_data2$Mean_GC)[[2]],
                  SemiLocalTrend = accuracy(pred_GC_semitrend_1$mean, new_data2$Mean_GC)[[2]])

GC_fit3_rmse <- data.frame(Commodity = "GC Fit 3",
                  DynamicRegression = accuracy(pred_GC_dynamic_3$mean, new_data3$Mean_GC)[[2]],
                  LinearTrend = accuracy(pred_GC_trend_3$mean, new_data3$Mean_GC)[[2]],
                  SemiLocalTrend = accuracy(pred_GC_semitrend_3$mean, new_data3$Mean_GC)[[2]])


combined_rmse <- rbind(CL_fit1_rmse, CL_fit2_rmse, CL_fit3_rmse, GC_fit1_rmse, GC_fit2_rmse, GC_fit3_rmse)

combined_rmse$Best_Model <- apply(combined_rmse[, -1], 1, function(x) {
  colnames(combined_rmse[, -1])[which.min(abs(x))]
})


print(combined_rmse)
##   Commodity DynamicRegression LinearTrend SemiLocalTrend  Best_Model
## 1  CL Fit 1          30.39734    6.255817      12.663403 LinearTrend
## 2  CL Fit 2          38.35451    3.689569      28.469102 LinearTrend
## 3  CL Fit 3          78.29768    2.932757       5.081877 LinearTrend
## 4  GC Fit 1         458.76358  249.223872     266.072375 LinearTrend
## 5  GC Fit 2         347.42962   11.774606     134.081005 LinearTrend
## 6  GC Fit 3        1531.99188  195.287713     204.353104 LinearTrend

QUESTIONS Local linear trend works best with everything except fit 3 of GC, which is best modeled using Semi Local Linear Trend state-specification. Local Linear Trend is best overall.

Lags are not useful in predicting CL and GC. All coefficients are 0 meaning no meaningful linear relationship between variables.